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Abstract 

For recursive circular filtering based on circular statistics, we introduce a general framework for 
estimation of a circular state based on different circular distributions, specifically the wrapped normal 
distribution and the von Mises distribution. We propose an estimation method for circular systems 
with nonlinear system and measurement functions. This is achieved by relying on efficient deterministic 
sampling techniques. Furthermore, we show how the calculations can be simplified in a variety of 
important special cases, such as systems with additive noise as well as identity system or measurement 
functions. We introduce several novel key components, particularly a distribution-free prediction 
algorithm, a new and superior formula for the multiplication of wrapped normal densities, and the 
ability to deal with non-additive system noise. All proposed methods are thoroughly evaluated and 
compared to several state-of-the-art solutions. 


1. Introduction 

Estimation of circular quantities is an omnipresent issue, be it the wind direction, the angle of 
a robotic revolute joint, the orientation of a turntable, or the direction a vehicle is facing. Circular 
estimation is not limited to applications involving angles, however, and can be applied to a variety 
of periodic phenomena. For example phase estimation is a common issue in signal processing, and 
tracking objects that periodically move along a certain trajectory is also of interest. 

Standard approaches to circular estimation are typically based on estimation techniques designed 
for linear scenarios that are tweaked to deal with some of the issues arising in the presence of circular 
quantities. However, modifying linear methods cannot only be tedious and error-prone, but also yields 
suboptimal results, because certain assumptions of these methods are violated. For example, solutions 
based on Kalman filters [1] or nonlinear versions thereof [2] fundamentally neglect the true topology 
of the underlying manifold and assume a Gaussian distribution, which is only defined on M”. In the 
linear case, the use of a Gaussian distribution is frequently justified by the central limit theorem. This 
justification no longer holds in a circular setting, as the Gaussian is not a limit distribution on the 
circle. 

In order to properly deal with circular estimation problems, we rely on circular statistics m, m , a 
subfield of statistics that deals with circular quantities. More broadly, the field of directional statistics 
[S] considers a variety of manifolds, such as the circle, the hypersphere, or the torus. Unlike standard 
approaches that assume linear state spaces, methods based on circular statistics correctly use the 
proper manifold and rely on probability distributions defined on this manifold. Gircular statistics 
has been applied in a variety of sciences, such as biology [1] , bioinformatics [B] , meteorology [7] , and 
geosciences |8]. 
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Table 1: Circular filters based on directional statistics. 


There has been some work on hltering algorithms based on circular statistics by Azmani et al. [2], 
which was further investigated by Stienne et al. m- Their work is based on the von Mises distribution 
and allows for recursive filtering of systems with a circular state space. However, it is limited to the 
identity with additive noise as the system equation and the measurement equation. The filter from [3] 
has been applied to phase estimation of GPS signals nn, m as well as map matching [TJ] . Markovic 
et al. have published a similar filter m based on the von Mises-Fisher distribution, a generalization 
of the von Mises distribution to the hypersphere. 

We have previously published a recursive filter based on the wrapped normal distribution allowing 
for a nonlinear system equation m- The paper extends this approach to make a nonlinear 
measurement update possible. Both papers rely on a deterministic sampling scheme, based on the 
first circular moment. This kind of sampling is reminiscent of the well-known unscented Kalman filter 
(UKF) We have extended this sampling scheme to the first two circular moments in so the 
proposed filters are, in a sense, circular versions of the UKF. 

The developed filters have been applied in the context of constrained tracking [IB], bearings-only 
sensor scheduling m. as well as circular model predictive control m- 

Furthermore, we proposed a recursive filter based on the circular Bingham distribution in m 
The Bingham distribution is closely related to the von Mises distribution, but has a bimodal density, 
which makes it interesting for axial estimation problems (i.e., problems with 180° symmetry). 

An overview of all of these filters and the considered distributions as well as system and measurement 
models is given in Table 

This paper summarizes and combines our results as well as extends the previous work by a number 
of additional contributions. First of all, we propose a general filtering framework that can be used in 
conjunction with a variety of system and measurement equations, different types of noise, and both 
the wrapped normal and the von Mises distributions. Our previous publications US], uni as well as 
the work by Azmani et al. jUj can be seen as special cases of the proposed framework. 

Furthermore, we introduce a new multiplication formula for wrapped normal distributions that 
outperforms the solution proposed in US). We generalize the prediction step from m to a purely 
moment-based solution that does not need to assume any kind of distribution. Compared to |Ii> 
we add the ability to deal with non-additive noise not only in the measurement update but also in 
the prediction step. Finally, we perform a thorough evaluation, where we compare the proposed 
techniques to several state-of-the-art approaches. 

This paper is structured as follows. First, we formulate the problem in Sec. Then, we introduce 
the necessary fundamentals from circular statistics in Sec.[^ Based on these fundamentals, we propose 
deterministic sampling schemes in Sec. and derive the operations on circular densities required for 
the circular filter in Sec. These results are used to introduce circular filtering algorithms in Sec. 
An evaluation of the proposed algorithms can be found in Sec. Finally, conclusions are given in 
Sec. m 

2. Problem Formulation 

In this section, we formulate the problems under consideration and summarize some standard 
approaches that have been used to address the issues associated with periodicity. 
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2.1. Circular Filtering 

Circular filtering considers estimation problems on the unit circle, which is commonly parameterized 
as the set of complex numbers with unit length {x G C : |a:| = 1}. To allow for a more convenient 
one-dimensional notation, we identify with the half-open interval [0,27r) C M, while keeping the 
topology of the circle. Together with the operation 

+ : X ^ , X + y := X -|-m y mod 27r, 

for all x,y € [0,27r) with standard addition -|-r on M, the circle forms an Abelian group. Because 
with the topology given above has a manifold structure, (5^, -|-) is a Lie group. 

We consider a system whose state Xk at time step /c is a value on the unit circle 5^. System and 
measurement models are assumed to be given. In this paper, we propose several methods to deal with 
different types of models. More complex models necessitate the use of more sophisticated algorithms 
and conversely, simpler models allow the use of computationally less expensive algorithms. 

2.1.1. System Model 

In this work, we consider a system model whose state evolves according to the general system 
equation 


Xk+l — 

with (nonlinear) system function ak ■ x W ^ and noise Wk € W stemming from some noise 
space W. Note that W is not necessarily S^, but may be an arbitrary set, for example the real-vector 
space M”, some manifold, or even a discrete set. An interesting and practically relevant special case is 
a (nonlinear) system with additive noise 

Xk+i = ak{xk) + Wk mod 27r 

with Ofc : 5*^ —)■ and Wk G S^. More particular, we also consider the special case, where is the 

identity, i.e.. 


Xk+i = Xk + Wk mod 27r. 


2.1.2. Measurement Model 

The system state cannot be observed directly, but may only be estimated based on measurements 
that are disturbed by noise. A general measurement function is given by 

^k ^fc(®A:Wfc) ) 

where Zk € Z is the measurement in the measurement space Z , hk ■ x V ^ Z is the measurement 
function and Vk € V is arbitrary measurement noise in a certain noise space V. Note that the 
measurement and noise space can be arbitrary sets in general. An interesting special case are 
measurement functions where the measurement noise is additive, i.e., 

Zk = hk{xk) -h Vk 

with measurement function hf^ : ^ Z and Vk £ Z. In this case, we require Z to have a group 

structure with -|- as the operation. Additionally, we consider the more specific case where hk is the 
identity, i.e.. 


Zk = Xk + Vk mod 27r 


with Zk,Vk G 5"^. 

Remark 1. We do not consider linear system models because linearity is a concept of vector spaces, 
not manifolds m- For this reason, there are no linear functions on the circle. 
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(a) Incorrect fusion. 

Figure 1: This figure illustrates that repositioning of the measurement is necessary to obtain satisfactory 
performance when using classical filters on circular problems. 



(b) Correct fusion. 
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2.2. Standard Approaches 

As circular estimation problems are wide-spread in a variety of applications, a number of standard 
approaches have been employed. We introduce three of the most common methods and explain their 
strengths and weaknesses. 

2.2.1. Gaussian-based approaches 

Gaussian-based methods (wrongly) assume a Gaussian distribution and use standard filtering 
techniques for Gaussians in conjunction with certain modihcations to allow their application to circular 
problems. 

One-Dimensional Methods. One common approach is the use of a standard Kalman filter [T] or, in 
case of nonlinear system or measurement functions, unscented Kalman filter (UKF) [2] with a scalar 
state Xk containing the angle 9k, i.e., Xk = Ok. However, two modifications are necessary before this 
approach can be used in practice. First, the estimate has to be forced to stay inside the interval [0,27r) 
by performing a modulo operation after every prediction and/or update step. 

Second, if the measurement space is periodic, the measurement needs to be repositioned to be 
closer to the state mean in certain cases. This problem occurs whenever the measurement and the 
current state mean are more than tt apart. In this case, the measurement needs to be moved by ±27rA: 
to an equivalent measurement that deviates at most tt from the state mean. An illustration of this 
problem is given in Fig. 

This type of filter is used as a comparison in m- When the uncertainty is small, this kind of 
approach works fairly well, but it tends to produce unsatisfactory results if the uncertainty is high. 

Two-Dimensional Methods. Another common approach is based on the Kalman filter or the UKF with 
two-dimensional state subject to a nonlinear constraint. More specifically, an angle 9k is represented 
by a state vector Xk = [cos(6lfc),sin(0)fc]^ and the constraint is ||xfc|| = 1 to enforce that is on the 
unit circle. In order to enforce this constraint, is projected to the unit circle after each prediction 
and/or update step. More sophisticated approaches increase the covariance to reflect the fact that the 
projection operation constitutes an increase in uncertainty m- 

This approach has been used in |18j . but did not perform as well as the filter based on circular 
statistics. One of the issues of this approach is the fact that the system and measurement model 
sometimes become more complicated when the angle 9k is transformed to a two-dimensional vector. 

2.2.2. Particle Filters 

Another method that can be applied is particle filtering [23] . Particle filters on nonlinear manifolds 
are fairly straightforward to implement because each particle can be treated separately. For the 
particle filter to work, the system function and the measurement likelihood both need to respect 
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(a) The wrapped normal distribution (red) is ob¬ 
tained by wrapping a Gaussian distribution (blue) 
around the circle. The parameters we used for this 
example are /r = 0, cr = 2. 



(b) The von Mises distribution (red) arises when 
restricting a two-dimensional Gaussian with /i = 
(cos /r, sin/i)^ and covariance k ■ 12 x 2 to the unit 
circle. The parameters for this example are /r = 

0, K = 1. 


Figure 2: Relation between WN and VM distributions and the Gaussian distribution. 


the underlying topology. The reweighting step as well as the commonly used sequential importance 
resampling (SIR) are independent of the underlying manifold and can be used in a circular setting as 
well. 

However, issues that are typically associated with particle filters arise. If the measurement 
likelihood function is very narrow, particle degeneration can occur, i.e., (almost) all particles have 
zero weight after the reweighting step. Furthermore, a large number of particles is required to obtain 
stable results. Even though these problems are less critical in a one-dimensional setting, there can 
still be issues if measurements with high certainty occur in areas with few particles. This can, for 
example, occur when information from sensors with very different degrees of accuracy is fused. It 
should also be noted that sampling from certain circular distributions can be somewhat involved and 
require the use of the Metropolis Hastings algorithm |24| or similar approaches, e.g., in case of the 
von Mises distribution. 


3. Circular Statistics 

Because of the drawbacks of the approaches discussed above, we propose a filtering scheme based 
on circular statistics [3], [3|- In the following, we introduce the required fundamentals from the field 
of circular statistics. 


3.1. Circular Distributions 

Circular statistics considers probability distributions defined on the unit circle. A variety of 
distributions has been proposed |1] . We give definitions of all distributions that are required for the 
proposed filtering scheme. 

Definition 1 (Wrapped Normal Distribution). The wrapped normal (WN) distribution is given by 
the probability density function (pdf) 


f{x]hW) 


1 

\/^a 


CO 

k=—oo 


(x 


fi -|- 2'Kk)‘^ 
" 2^2 


with X G S^, location parameter pL G , and concentration parameter cr > 0. 
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The WN distribution is obtained by wrapping a one-dimensional Gaussian distribution around 
the unit circle and adding up all probability mass that is wrapped to the same point (see Fig. 2a). 
It appears as a limit distribution on the circle m in the following sense. A summation scheme of 
random variables that converges to the Gaussian distribution in the linear case, will converge to the 
WN distribution if taken modulo 27r. Because of its close relation to the Gaussian distribution, the 
WN distribution inherits a variety of its properties, for example its normalization constant as well as 
the formula for the convolution of densities. Even though there is an infinite sum involved, evaluation 
of the pdf of a WN distribution can be performed efficiently, because only few summands need to be 
considered 


Definition 2 (Von Mises Distribution). The von Mises (VM) distribution is given by the pdf 

f{x;n,K) = ] exp(Kcos(x - p)) 

with X £ S^, location parameter p, £ , and concentration parameter n > 0. Io{-) is the modified 

Bessel function of order 0. 

The VM distribution, sometimes also referred to as the circular normal (GN) distribution, is 
obtained by restricting a two-dimensional Gaussian with mean ||/r|| = 1 and covariance k ■ 12x2 (where 
I 2 x 2 is the identity matrix) to the unit circle and reparameterizing to [0,27r) as can be seen in Fig. 

For this reason, it also inherits some of the properties of the Gaussian distribution; most importantly, 
it is closed under Bayesian inference. The VM distribution has been used as a foundation for a circular 
filter by Azmani et al. |U] . It is closely related to the Bingham distribution and conversion between the 
two is effectively a matter of shrinking the the interval [0,27r) by a factor of two two, i.e., /(2x; p, k) 
is a Bingham distribution m p- 4]. 

Definition 3 (Wrapped Dirac Mixture Distribution). The wrapped Dirac mixture (WD) distribution 
with L components is given by 


L 

/(x; 7 i,... ,7 l,/3i, ... ,/3l) = ^7j(5(x - fifi 


with Dirac delta distribution 6{-), Dirac positions fii,..., /3 l £ S^, and weights 71 ,... ,7l >0 where 
Et=i7i = l- 

Unlike the WN and VM distributions, the WD distribution is a discrete probability distribution 
on a continous domain. It can be obtained by wrapping a Dirac mixture in M around the unit circle. 
WD distributions can be used as discrete approximations of continuous distributions with a finite set 
of samples. 

In this paper, we use the following notation. We denote the density function of a WN distribution 
with parameters p and a with WAf{p, cr). If a random variable x is distributed according to this WN 
distribution, we write x ~ WM{p, a). The terms VM.{p, k) and >VT’( 7 i,..., 7 l, fii,..., fdi) are used 
analogously to describe the density functions of VM and WD distributions with parameters p, k and 
71 ,..., 7 l, /3i,..., /3 l, respectively. 


3.2. Circular Moments 

In circular statistics, there is a concept called circular (or trigonometric) moment. 

Definition 4 (Gircular Moments). For a random variable x ~ f{x) defined on the circle, its n-th 
circular moment is given by 

mn = E(exp(zx)"') = E(exp(znx)) 
exp(mx) • f{x) dx G C 

with imaginary unit i. 
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The n-th moment is a complex number and, hence, has two degrees of freedom, the real and 
the imaginary part. For this reason, the first moment includes information about the circular mean 
arg 772-1 = atan2(Im mi , Re ttt-i) as well as the concentration | mi\ = \J (Re?72i)^ + (Im 7721 )^ of the 
distributiorQ similar to the hrst two real-valued moments. R can be shown that WN and VM 
distributions are uniquely defined by their first circular moment [H]. 

Remark 2. Any continuous and piecewise continuously differentiable 27r-periodic function / : M —)■ M 
can be written as a Fourier series 


f{x) = ^ Cfc ex.Y>{ikx) 


k=—oo 


r27T 


where Ck = 


27r 


f{x)exp{—ikx)dx 


ff f{') Is the pdf of a circular probability distribution, the Fourier coefficients are closely related to the 
circular moments according to c^ = ^r77_fc. 

Lemma 1 (Moments for WN, VM, and WD Distributions). For WN, VM, and WD distributions 
with given parameters, the n-th circular moment can be caleulated according to 


WN 

772„ = exp ( 272/i 




m]^^ = exp(m/x) 


2 

^\n \(^) 

Io{k) 


777, 


WD 


= ^7jexp(m/3j) . 


A proof is given in [S] . 


3.3. Circular Moment Matehing 

As both WN and VM distributions are uniquely dehned by their first moment, it is possible to 
convert between them by matching the first circular moment. 


Lemma 2 (Circular Moment Matching). We define A(x) 


as given in 


1. For a given first moment mi, the WN distribution with this first moment has the density 
WN ^atan2(Im7721, Re 7771), -\/—21og (|7?7i|)^ . 

2. For a given first moment mi, the VM distribution with this first moment has the density 
VM (atan2(Im 777i,Re 772i), A“^(|777 i|)). 

3. For a given VM distribution with density VM{p, k), the WN distribution with identical first 
moment has the density WAf ( p,, 



4. For a given WN distribution with density WJ\f{yL,a), the VM distribution with identical first 


moment has the density VM (^pi,A ^ ^exp - 


The proof is given in m- Calculation of the function A ^(O is somewhat tricky. In m , we use 
the algorithm by Amos |28j^ to calculate A(-) and MATLAB’s fsolve to invert this function. Stienne 
et al. have proposed closed-form approximations, which can be calculated very easily but have a large 
approximation error m, [in]. A more detailed discussion on approximations of A ^(•) can be found 
in |29| and |30l Sec. 2.3]. 


^The term 1 — |mi| is sometimes referred to as circular variance. 
^Pseudocode of this algorithm is given in m- 
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4. Deterministic Sampling 


In order to propagate continous probability densities throngh nonlinear fnnctions, it is a common 
technique to nse discrete sample-based approximations of the continons densities. A set of samples can 
easily be propagated by applying the nonlinear function to each sample individually. This approach 
can be used for both the prediction and the measurement update step. 

We distinguish between deterministic and nondeterministic sampling. Nondeterministic sampling 
relies on a randomized algorithm to stochastically obtain samples of a density. Typical examples 
include the samplers used by the particle filter [23] or the Gaussian particle filter m- Deterministic 
sampling selects samples in a deterministic way, for example to fit certain moments (the sampler 
used by the UKF, El), or to optimally approximate the shape of the density (published in j32j . this 
sampler is used by the S^KF, [ 23 ] )• Deterministic sampling schemes have the advantage of requiring a 
significantly smaller number of samples, which is why we will focus on this type of solution. 

A naive solution for approximating a WN density may be the application of a deterministic 
sampling scheme for the Gaussian distribution (such as the samplers used in [2] , [ 33 ] ) and subsequently 
wrapping the samples. Even though this technique is valid for stochastic samples, it does not provide 
satisfactory results for deterministic samples. In extreme cases, wrapping can cause different samples 
to be wrapped to the same point, grossly misrepresenting the original density. This problem is 
illustrated in Fig. In the case of UKF samples (Fig. 3a), it can be seen that for a ~ 2.5, one 
sample is placed at /r and two samples are placed on the opposite side of the circle, i.e., the mode 
of the approximation is opposite to the true mode. Furthermore, for cj ~ 5 all three UKF samples 
are wrapped to the same position, i.e., the sample-based approximation degenerates to a distribution 
with single Dirac component even though the true distribution is nearly uniform. Similar issues arise 
in the case of S^KF samples (see Fig. 3b I. 


4 . 1 . Analytic Solutions 

First of all, we consider analytic solutions to obtain deterministic samples. These solutions are 
based on circular moment-matching and only provide a small, fixed number of Dirac delta components, 
but are extremely fast to calculate, making them a good choice for real-time applications. 

In |15] . we presented a method to obtain a WD approximation with three equally weighted 
components, which is based on matching the first circular moment (see Algorithm [^. We further 
extended this scheme to obtain a WD with five components by matching the first as well as second 
circular moment (see Algorithm , which, as we proved in necessitates the use of different 
weights. Both approaches can approximate arbitrary symmetric circular densities with given circular 
moments. 


Algorithm 1: Deterministic approximation with L = 3 components. 

Input: first circular moment mi 
Output: WT>(7 i,...,73,^i,...,^3) 

/* extract */ 

fj, •(— atan2(Immi,Remi); 

/* obtain Dirac positions */ 

a <— arccos (||mi| — 2 ); 
j3i yi — a mod 27r; 

/32 •(— /r mod 27r; 
k + ot mod 27r; 

/* obtain weights */ 

71,72,73 ^ 5; 


8 







Algorithm 2: Deterministic approximation with L = 5 components. 

Input: first circular moment mi, second circular moment m 2 , 
parameter A e [0,1] with default A = 0.5 
Output: >VP( 7 i,... , 75 ,^ 1 ,... ,/ 35 ) 

/* extract ^ */ 

fj, •(— atan2(Immi,Remi); 
mi |mi|; 
m 2 |m 2 |; 

/* obtain weights */ 

7 ™“ (4mf — 4mi — m2 + l)/(4mi — m 2 — 3); 

7 max ^ (2mf — m 2 — l)/(4mi — m 2 — 3); 

75 ^ 75”“ + A( 75 ”“ - 75”“); 

71,72,73,74 ^ (1 - 75)/4; 

/* obtain Dirac positions */ 

Cl ^ i^(mi - 75 ); 

C 2 ^ l^("l2-75) + l; 

X 2 (2ci + 

Xl Cl - X2] 

(j)! ^ arccos(xi); 

4>2 ^ arccos(x 2 ); 

(/3i,... ,/35) ^ /I + -4>2,+<p2,0) mod 27r; 



0 5 


(a) Analytic solution from Algorithm and 
wrapped UKF [5] samples. 



0 5 


(b) Analytic solution from Algorithm with A 
0.8 and wrapped S^KF samples. 


Figure 3: Proposed approaches for generating samples of WN distributions with a different concentra¬ 
tion parameters a compared to the naive approach of wrapping samples of a Gaussian with identical 
a. It can be seen that the UKF and S^KF samples are eventually wrapped to the same location, 
which produces an extremely poor approximation. 
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Figure 4: Example of the deterministic sampling of a VM and WN distribution with equal first circular 
moment. From top to bottom: a) original densities, b) result of Algorithmic c) result of Algorithmic 
d) approach based on VM kernels from [31] for 10 components, e) quantization approach from |3S] for 
10 components. Note that the result of Algorithm |C is identical for both densities because only the 
first circular moment is considered. 


4-2. Optimization-based Solutions 

If a larger number of samples is desired and there are more degrees of freedom in the samples 
than constraints (such as preservation of circular moments), optimization-based solutions can be used. 
The number of samples can be adjusted by the user and an optimal approximation is derived by 
minimizing a distance measure. 

In order to simultaneously calculate optimal locations and weights for the samples, a systematic 
approach based on VM kernels has been proposed in [31]. For a WD mixture, an induced VM mixture 
is compared to the true distribution with a quadratic integral distance. A specific kernel width is 
considered for each component, which depends on the weight of the component and the value of the 
true distribution at the location of the component. Both the weights and the locations of a fixed even 
number of WD components are optimized to obtain an optimal symmetric approximation. Constraints 
in the optimization algorithms are used to maintain a predefined number of circular moments. This 
approach results in well-distributed Dirac mixtures that fulhll the moment constraints. 

A quantization approach is discussed in [35] . It is based on computing optimal Voronoi quantizers. 
In this approach, optimality refers to minimal quadratic distortion. The resulting Voronoi quantizer 
gives rise to a circular discrete probability distribution on a continuous domain that approximates the 
original continuous distribution. Use of this approximation is particularly beneficial in the prediction 
step of stochastic filters, because an error bound for propagation through a non-trivial system function 
can be obtained without actually knowing that function. It is sufficient to require it to be Lipschitz 
and to know an upper bound for the Lipschitz constant. Furthermore, circular moment constraints 
can be introduced in the optimization procedure of the quantization approach. 

Examples from all discussed methods for deterministic sampling are depicted in Eig. |^ 

5. Operations on Densities 

In order to derive a circular filtering algorithm, we need to be able to perform certain operations 
on the involved probability densities. 

5.1. Shifting and Mirroring 

For a given density /(x), we want to obtain the density f{c — x) for a constant c G S^. This 
operation is necessary in certain cases of the update step. We can split this operation into two steps: 
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mirroring to obtain /(—x), and subsequent shifting by c to obtain /(c+ (—x)). Mirroring WM{n,a) 
and VA^(/U,k) yields WM{2 'k — fi,a) and VA^(27r — fi,K) because the distributions are symmetric 
around their mean. Shifting and by c yields 

WM{fi — c mod 27r, cr) and — c mod 27r, k) , 

so the combined operation results in 

WAA((27r — fj.) — c mod 27r, a) 
and V2Vl((27r — fi) — c mod 27r, k) . 


5.2. Circular Convolution 

Given two independent circular random variables xi r\j flixi),X 2 / 2 (x 2 ), the sum xi + X 2 is 
distributed according to 


(/i * f2){x) 


2tt 


flit)f2ix 


t) dt , 


Jo 

where * denotes the convolution. This operation is necessary in the prediction step to incorporate 
additive noise. 

WN distributions are closed under convolution and the new pdf can be obtained just as in the 
Gaussian case m, i-e., /r = /ii + //2 mod 27r, = af + cr^. VM distributions are not closed under 

convolution. For this reason, Azmani et al. |9] used the approximation from [5], which is given by 
/X = /ri + /U 2 ,k = A~^{A{k,i), A{k 2 )). The function A(-) is the same as defined in Lemma This 
approximation can be derived from an intermediate WN representation m- A similar approximation 
has been used by Markovic et al. for the von Mises-Fisher case jl4l (7)]. 

In this paper, we present a more general result that calculates the convolution based on circular 
moments. 


Lemma 3 (Moments After Addition of Random Variables). Assume independent random variables 
Xl ~ /l,X2 f 2 defined on the circle. For the sum x = xi + X2, it holds 

E(exp(inx)) = E(exp(inxi))E(exp(inx 2 )) . 


Proof. 


n27Z 

mn =E(exp(mx)) = / exp(inx)/(x) dx 

Jo 

/*27r /*27r 

= / / exp(m(x))/i(y)/2(x - y)dydx 

Jo Jo 

r27T p27r 

= / exp(m(xi+ X2))/i(xi)/2(x2)dxidx2 

Jo Jo 

p27T p2n 

= / exp(inxi)/i(xi) dxi / exp(inx2)/2(a;2) da:2 

Jo Jo 

=E(exp(mxi))E(exp(inx 2 )) 


□ 


If moment matching of the first circular moment is used to fit a WN or VM to the density that 
results from convolution, the solutions for WN and VM distributions from [9] and m arise as special 
cases of Lemma [3l 

Remark 3. In fact, Lemma\^ allows us to calculate the convolution purely moment-based. For this 
reason, we do not need to assume any particular distribution, but can just calculate the moments of 
the convolved density based on the produets of original moments^ 

®In general, for example in the case of mixture densities, the complexity of the density increases with each successive 
convolution, so considering only a finite number of moments constitutes an approximation. 
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0.7 



Figure 5: Multiplication of two WN densities with parameters //i = 2 ,(Ti = 0.7, and 
fi 2 = 4.95,172 = 1-3. The true product and the results of both proposed approximation methods 
(VM and moment-based) are depicted. Note that the true product is not a WN density. 


5.3. Multiplication 

Multiplication of pdfs is an important operation for filtering algorithms, because it is required for 
Bayesian inference. In general, the product of two pdfs is not normalized and thus, not a pdf. For 
this reason, we consider the renormalized product, which is a valid pdf. 

5.3.1. VM 

Von Mises densities are closed under multiplication |^. It holds that VA^(/xi, ki) ■ ^ 2 ) oc 

k), where 


/i = atan2(Immi,Remi), K=\mi\, 
with mi = Ki exp{ifii) + K 2 exp(i^ 2 ) ■ 


5.3.2. WN 

WN densities are not closed under multiplication. In the following, we consider two different 
methods to approximate the density of the product with a WN density. 

WN via VM. WN densities are not closed under multiplication. In we proposed a method to use 
the VM distribution in order to approximate the product of two WN densities. More specifically, we 
convert the WN densities to VM densities using Lemma multiply according to the VM multiplication 
formula, and convert back to a WN distribution by applying Lemma again. This method has the 
disadvantage that, in general, the first circular moment of the resulting WN does not match the first 
circular moment of the true product. An example can be seen in Fig. 

WN via Moment Matching. In this paper, we present a new method for approximating the product of 
WN distributions. This method is based on directly approximating the true posterior moments. 

Theorem 1. The first circular moment of WAf{fJ,i, ai) ■ WM[pL 2 ,(J 2 ) after renormalization is given 
by 


mi 


00 277 

E w{j,k) f exp(ix)Af(x; /i(j, k), a(j, k))dx 

j\k=—oo 0 

CXD 277 

E w{j,k) f M{x-,p.{j,k),a{j,k))dx 

j,k=—OQ 0 
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where M{x] fx, a) is a one-dimensional Gaussian density with mean fx and standard deviation a, and 


k) 


w{j, k) 


{jjLi + 27rj)fT2 + (/i2 + 2TTk)aj 


o-f+ 0-2 


2^2 


at a. 


1^2 


o-f+ cj| 


exp 


1 ((/Ji+27rj)-(^2+27rfc))^ 

2 crj^+o-l 


y/27r{a‘( + a‘i) 


Proof. The true renormalized product is given by f{x) = c ■ f{x]fxi,ai) ■ f{x]fX 2 ,a 2 ), where c 
renormalizes the product, i.e., 


c = 



f{x-,fxi,ai) ■ f{x-,fX 2 ,a 2 )dx 


-1 


We calculate 


r2n 


mi=c- / exp{ix) ■ f{x] fxi,ai) ■ f{x-, iX 2 ,a 2 )dx 


JO 

f' 2 'K 

=c ■ / ex 

Jo 


P(^a;) • ^ M{x] fxi + 27rj, cji) 
j=-oo 

CO 

AA(x; //2 + 27rA:, CT 2 ) dx 


k=—oo 


CO CO p27V 

=c • W W / exp(ix) • AA(x; ;Ui + 27 rj, cJi) 

■ JO 
J = — OQ k= — oo 

■ Af{x] fX2 + 2 'Kk, (T2) dx 

CO CO p27V 

=C- / exp(ix) • u;(j, /c) 

j=—oo k=—oo ® 

• •AA(x;^(i, /c),a(j,A:))dx 

CO CO p27V 

=c • V V -wlyj, k) ■ / exp(ix) 

- 7 - 

j=—oo k=—oo 

■Af{x;ix{j,k),a{j,k))dx , 


where we use the dominated convergence theorem to interchange summation and integration. We use 
the abbreviations, 


P(j, k) 
o-U, k) 


w{j, k) 


(^1 + 2 T:j)al + {fX 2 + 2 'nk)a\ 

al + al 


2^2 


at a. 


1^2 


af-ha^ 


exp 


1 ((;ji+27rj)-(;j2+27rfc))^ 

2 o-f+o-l 


v/MofTfrf) 


based on the multiplication formula for Gaussian densities given in j36[ 8.1.8]. 
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To obtain the renormalization factor c we use a similar derivation 

/■27r 


c-' = 


■ f{x]^i2,a2)dx 


10 


p27T ^ 

■- / ^ AA(x;/xi + 27rj,c7i) 

0 j=-oo 

OO 

■ ^ AA(x;/i2 + 27rfc, 1 T 2 ) dx 


k=—oo 
00 00 

E 

j=—oo k=—oo ' 


r27T 


M{x]fii +27rj, cJi) 


• AA(x; fi2 + 27rA:, 1 T 2 ) dx 


CO CO 


: Y '^U,k) 

j=—oo k=—oo 
p27V 

■ / ■^{x-,fj.{j,k),a{j,k))dx 


□ 


The involved integrals can be reduced to evaluations of the complex error function erf [261 7.1]. 
This yields 

^277 

/ exp(ix) • A7(x; n{j,k),a{j,k)) dx 

Jo 

= \ exp k) - 


erf 
— erf 


Kj,k) + icf{j,kf 
V2a{j, k) 
fi{j,k) -27r+ icr{j,kf 
k) 


and 


r-27r 


J^{x]fi{j,k),a{j, k))dx 


= i lerf f 

2 V \a{j,k)y/2 


— erf 


k) - 27r 
a{j,k)y/2 


There are efficient implementations of the complex error function by means of the related Faddeeva 
function m- Furthermore, the infinite sums can be truncated to just a few summands without a 
significant loss in accuracy. For example, the multiplication in Fig. [^requires 5x5 summands for 
an error smaller than the accuracy of the IEEE 754 64 bit double data type [SH]. Consequently, the 
proposed method allows for efficient calculation of the approximate multiplication of WN densities. 

6. Circular Filtering 

Based on the results in the previous section, we derive circular filtering algorithms for the scenarios 


described in Sec. 2.1 All proposed algorithms follow the recursive filtering concept and consist of 
prediction and measurement update steps. We formulate the necessary steps without requiring a 
particular density whenever possible such that most methods can be directly applied to WN as well 
as VM distributions, and might even be generalized to other continous circular distributions. An 
overview of all considered prediction and measurement update algorithms is given in Table 
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system model 



measurement 

model 

method 

identity 

additive noise 

arbitrary noise 

identity 

additive noise 

arbitrary noise 

WN 

VM 

moment-based 

1151 (special case) 

contribution 

[m 

contribution 

contribution 

contribution 

contribution 

contribution 

i 

\M 

contribution 

EH 

contribution 


Table 2: Prediction and measurement update algorithms. Entries marked with contribution are 
contributions of this paper. 


6.1. Prediction 

The prediction step is used to propagate the estimate through time. 

6.1.1. Identity System Model 

The transition density is given according to 

I'2'k 

f{Xk+l\Xk) = / f{Xk+l,Wk\Xk)<IWk 


r-27r 


= / f{Xk+l\Xk,Wk)f'“{Wk)&.Wk 

Jo 

l' 2 n 

= / 5{xk+i-{xk + Wk))f'^{wk)dwk 

Jo 

= /’^(Xfc+1 - Xk) , 

where /"'(•) is the density of the system noise. For the predicted density, according to the Chapman- 
Kolmogorov equation we obtain 


r 2 TT 


fP{xk+i) = / f{xk+i\xk)nxk)dxk . 


( 1 ) 


In the special case of an identity system model, this yields 


r2n 


fP{xk+i)= / f^{xk+i - Xk)f{xk)dxk 

Jo 

= {r*n{xk+i), 


where * denotes convolution as defined in Sec. 5.2, For the VM distribution, this system model has 
been considered in [2] . If a WN distribution is assumed, Q is a special case of m where we omit the 
propagation through the nonlinear function. 

6.1.2. Nonlinear System Model with Additive Noise 

Similar to the previous case, the transition density is given by 

p2'K 

f{xk+i\xk) = / d{xk+i - {ak{xk) + Wk))f'^{wk) dwk ■ 

Jo 

We approximate the prior density f^{xk) ~ J2j=ilj^ixk — fJj) using, for example. Algorithmor 
Algorithm Then, the prediction density can be approximated according to 


Fixk+i) 


r 2 TT 


f{xk+i\xk)fixk) dxk 
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/ /*27r 


'0 


5{xk+i - {akixk) + Wk))r{wk) ddwk 
^7j(5(xfc -/3j)j dxfc 

^27r 

/ (^(xfc+i - (afe(xfc) + Wfe)) 
,=i -^0 


^i=i 

r-27r 


5{xk - I3j)dxk dwk 


r-27r 


f'“{wk) ^ 7j(^(xfc+i - (afc(^j) + Wk)) dwk 


i=i 


^fi^k + l-Wk) 


= (/"' * /)(a^fe+i) , 


where / is obtained by moment matching of the WD Yli= i 7j<5(x — {ak{(3j)) according to Lemma |2 
The convolution can be calculated as described in Sec. 15.21 


6.1.3. Nonlinear System Model with Arbitrary Noise 

In this paper, we extend the previous results to deal with arbitrary noise in the prediction step. 
For arbitrary noise, the transition density is given by 

p2n 

f{xk+i\xk)= / 5{xk+i-ak{xk,Wk))f'^{wk)dwk . 

Jo 

We approximate the prior density f^{xk) ~ Yl!j=i — fJj) as well as the noise density /"'(rcfc) 

Yli=ilT— fJf)- It should be noted that the noise is not necessarily a circular quantity and 
different approximation techniques may be required. If IF = M"', the techniques presented in [35] may 
be applied. Then, the prediction density can be approximated according to 


F{xk+i) 

r 2 TT 


f{xk+i\xk)r{xk)dxk 


r 2 TT 


d{xk+i - ak{xk, Wk))f'^{wk) dwkfixk) dxk 


/o JW 

r2iT 


6 {xk+i- ak{xk,Wk)) PD 


-. 1=1 


JO JW 
dwk f{xk) dxk 
1 - 2 ^ 

= y^7P / ^(^k+i - ak(xk, PP))f{xk) dxk 

1=1 

T rk 

^ l' 2 n 

/ ^i^k+i - ak{xk, PT)) 

1=1 

'^Ijdixk-Pj^ dxk 

L L™ 

EE ljlT^{xk+i - ak{Pj,PP)) 

j=i 1=1 
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A continuous density can be fitted to this result by circular moment matching. The algorithm is given 
in Algorithmic It is worth noting that this algorithm can be executed based purely on the circular 
moments of f{xk) and f^{wk) if the deterministic sampling scheme only depends on these moments. 
In that case, we do not necessarily need to fit a distribution f{xk+i) to the resulting circular moments, 
but can just store the estimate by retaining those circular moments. 


Algorithm 3: Prediction with arbitrary noise. 

Input: prior density f{xk), system noise density f'^{wk), system function ak{-,-) 

Output: predicted density f{xk+i) 

/* sample prior density and noise density */ 

(71,... ,7l,/3i, ... ,/ 3 l) ^ sampleDeterm(/(a:A:)); 

( 7 ^, • • • ^ sampleDeterm(/“'(u;fc)); 

/* obtain Cartesian product and propagate */ 

for j ^ 1 to L do 
for / 1 to L"' do 

^j+L(i-i) ^ • 7/"; 

P^+L{1-1) ^ Pf’)', 

end 

end 

/* obtain posterior density */ 

f{xk+i) ^ momentMatching( 7 ^,... , 7 £.^„, /3f,..., 


6.2. Measurement Update 

The measurement update step fuses the current estimate with a measurement that was obtained 
according to the measurement equation. 

6.2.1. Identity Measurement Model 

In the case of the identity measurement model and additive noise, the measurement likelihood is 
given by 


f{Zk\Xk) = f{zk - Xk) . 

For the posterior density, application of Bayes’ theorem yields 

f{Xk\Zk) = 0 ^ f{Zk\Xk)f{Xk) , 

where f{xk) is the prior density. Thus, we obtain the posterior density 

f{Xk\Zk) OC f{Zk - Xk)f{Xk) , 


as the product of the prior density and f^izk — Xk), which can be obtained as described in Sec. 5.1 


The multiplication depends on the assumed probability density and can be performed using the 


multiplication formulas given in Sec. 5.3 For the VM case, this is equivalent to the measurement 
update from [9] and for the WN case, this is equivalent to the measurement update from m- 


6 . 2 . 2 . Nonlinear Model with Additive Noise 

For a nonlinear measurement function with additive noise, the measurement likelihood is calculated 
according to f{zk\xk) = P{zk — hk{xk)), as given in (TU]. The remainder of the measurement update 
step is identical to the case of arbitrary noise as described in the following section. 
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6.2.3. Nonlinear Model with Arbitrary Noise 

For a nonlinear update with arbitrary noise, we assume that the likelihood is given. The key idea 
is to approximate the prior density f{xk) with a WD mixture and reweigh the components according 
to the likelihood f{zk\xk). However, this can lead to degenerate solutions, i.e., most or all weights are 
close to zero, if the likelihood function is narrow. 

We have shown in that a progressive solution as introduced in [32] can be used to avoid this 
issue. For this purpose, we formulate the likelihood as a product of likelihoods 

fUklXk) = fiZk\Xk)^^ • ... • fiZklXk)^’’ , 

where Ai,..., > 0 and This decomposition of the likelihood allows us to perform the 

measurement update step gradually by performing s partial update steps. Each update step is small 
enough to prevent degeneration and we obtain a new sample set after each step, to ensure that the 
differences between the sample weights stay small. In order to determine Ai,..., A^ and s, we require 
that the quotient between the smallest weight 7min and the largest weight 7max after reweighing is 
not below a certain threshold R G (0,1), i.e., > ji Using the conservative bounds 

7min > . min ( 7 ”) ' . min (f{zk\/3^)) , 

3=1,...,L 3=1,■■■,L 

7 max < . niax ( 7 ^) • niax (/(zfc |/3”)) , 

3=1,...,L 3=1, 


this leads to the condition 


R- ^ 

3 

log J 



max ( 7 ") 
min ( 7 ") 


where WI 1 ( 7 ”,..., 7 ^, ,0]^,..., /3£) is the deterministic approximation at n-th progression step|^ The 
progression continues until A^ = 1. This method can be applied in conjunction with WN as well 
as VM distributions (see Algorithm]^. 


7. Evaluation 

7.1. Propagation Aeeuraey 

In order to evaluate the deterministic sampling as introduced in Sec. we investigate the accuracy 
when performing propagation through the nonlinear function g : ^ 

g{x) = X + c X]R sin(x) mod 2tt , 

where c G [ 0 , 1 ) is a parameter controlling the strength of the nonlinearity and Xr refers to multipli¬ 
cation in the field of real numbers M. Furthermore, we consider the density WAA(0, a) that we want to 
propagate through g(-). For this purpose, we sample >VAA(0,cr) deterministically using the methods 
described in Sec. [^and obtain >VP( 7 i,..., 7 ^, Pi,..., Pl)- Then, we apply g{-) componentwise, which 
yields WVI^i,.. .,-iL,g{Pi), ■ ■ • , 5 (^l))- 
The true posterior is given by 

^ f {x); a) 
g'{x) 


^Compared to m, we extend the progressive scheme to handle discrete approximations with non-equally weighted 
components. 
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Algorithm 4: Progressive measurement update for arbitrary noise. 

Input: measurement likelihood predicted density f^{xk) as WN or VM, threshold 

parameter R 

Output: estimated fensity as WN or VM 

s 0; 

f{xk) ^ Fixk)] 

while Y^n=i An < 1 do 
s i — S H" 15 

/* deterministic sampling from WN or VM density (Sec. */ 

( 71 ,... ,7l,/3i, ... ,/3l) ^ sampleDeterm(/(xfc)); 

/* calculate size of progression step */ 

R\ \ 


A. 


mm 


1 - a. 


log R 


log 


max {7,') 
min (7 ) 




/* execute progression step 
for j ^ 1 to L do 

end 

/* use moment matching to obtain WN or VM density 

/ ^momentMatching( 7 i,..., 7 ^, /3i,..., ^l); 


*/ 


*/ 


end 


and can only be calculated numerically. We evaluate the first and the second circular moment 
mWD^ i = 1,2 of the resulting WD distribution and compare to the hrst and the second circular 
moment i = 1,2 of the true posterior, which is obtained by numerical integratiorj^ The 

considered error measure is given by i = 1,2, where | • | is the Euclidean norm in the 

complex plane. Additionally, we ht a WN density to the posterior WD by circular moment matching 
and numerically calculate the Kullback-Leibler divergence 


Dkl 


l•2n / ftrue/^\ 


( 2 ) 


between the true posterior and the htted WN. The Kullback-Leibler divergence is an information 
theoretic measure to quantify the information loss when approximating with This concept 

is illustrated in Fig. 

The results for different values of a are depicted in Fig. We compare several samplers, the 
analytic methods with L = 3 components (Algorithm and L = 5 components (Algorithm 
parameter A = 0.5) from Sec. 4.1 as well as the quantization approach discussed in Sec. 


4.2 It can 


be seen that the analytic solution with L = 5 components is signihcantly better than the solution 
with L = 3 components. The quantization-based solution is computationally quite demanding but 
gives almost optimal results. However, the analytic solution with L = 5 components has comparable 
performance in terms of the Kullback-Leibler divergence even though the posterior moments are not 
calculated as accurately. 

1.2. Moment-Based WN Multiplieation 

In this evaluation, we compare the two methods for WN multiplication given in Sec. |5.3.2| and 
For two WN densities and we calculate the true product 

(Ti) • yVAf{p. 2 W 2 ) and compare it to the WN approximation 


Sec. 5.3.2 

jtrue ^ 


^Numerical integration produces very accurate results in this case, but is too slow for use in practical filtering 
applications. 
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prior g(x) = X + csm{x) mod 27r posterior 

Figure 6: Propagation of a WN distribution with parameters g, = 0.1, cr = 1 through a nonlinear 
function g by means of a deterministic WD approximation with five components. In this example, we 
use c = 0.7. 


In order to determine the approximation quality, we compute the Kullback-Leibler divergence 
Furthermore, we consider the distance 


^fitted 


r- 27 r 


f^ftrue^x) - dx . 


The results for different values of cri, fj, 2 , and a 2 are depicted in Fig. |^and Fig. We keep gi fixed 
because only the difference between g 2 and ni affects the result. Multiplication is commutative, so we 
consider different sets of values for ui and a 2 to avoid redundant plots. 


As can be seen, the moment-based approach derived in Sec. 5.3.2 signihcantly outperforms the 


approach from Sec. 5.3.2 in almost all cases according to both distance measures. The new approach 
is particularly superior for small uncertainties. 


7.3. Filtering 


Scenario System Function Measurement Noise C’' 


s 

m 

1 

s-non-additive 
m-non-additive 
1-non-additive 


(3) 

(3) 

(4) 
(4) 
(41 


0.01 

I2x2 

0.1 

I2x2 

3 

I2x2 

0.01 

I2x2 

0.1 

I2x2 

3 

I2x2 


Table 3: Evaluation scenarios. 


In order to evaluate the proposed filtering algorithms, we simulated several scenarios. First of all, 
we distinguish between models with additive and with a more complex noise structure. In the case of 
additive noise, we consider the system function 

Xk+i =Xk + Cl Xr sin(xfe) +C 2 + Wk , (3) 

with two parameters ci = O.I,C2 = 0.15, noise Wk ~ WA7(0,0.2), and Xr is multiplication in the field 
of real numbers M. Intuitively, ci determines the degree of nonlinearity and C 2 is a constant angular 
velocity that is added at each time step. For the case of arbitrary noise, the system function is given 
by 


Xk+l =Xk + Ci Xr sin(xfc -f Wk) + C2 , 


(4) 
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a = 0.5 


cr = 1.0 
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Figure 7: Propagation of WA7(0 ,(t) through a nonlinear function with nonlinearity parameter c. 


with the same ci, C 2 , and as above. In both cases, the nonlinear measurement function is given by 

h = [cos(xfc),sin(xfc)]^ + :Ufc £ 

with additive noise ~ A7(0 ,77 • 12 x 2 ), V £ {3,0.1,0.01}. An overview of all considered scenarios is 
given in Table 

In the scenarios with additive system noise, we compare the proposed filter to all standard 
approaches described in Sec. |2.2[ a UKF with ID state vector, a UKF with 2D state vector and particle 
filters with 10 and 100 particles. In order to handle non-additive noise with a UKF, typically state 
augmentation is used, which is not applicable to arbitrary noise. For this reason, we only compare the 
proposed approach to the particle filters in the non-additive noise case. 

The initial estimate is given by xq ~ WA7(0,1), whereas the true initial state is on the opposite 
side of the circle = vr, i.e., the initial estimate is poor, which is difficult to handle for noncircular 
filters. 

For the circular filtering algorithm, we use the deterministic sampling method given in Algorithm]^ 
with parameter A = 0.5. The progression threshold is chosen as R = 0.2. 

In order to evaluate the performance of different filters, we consider a specific error measure that 
takes periodicity into account. The angular error is defined as the shortest distance on the circle 

d : X ^ [0, vr] , d{a, h) = min(|a — 6|, 27r — |a — b\) . 

This leads to an angular version of the commonly used root mean square error (RMSE) 


1 



max 


^max 

^ d{xk, 

k=l 


„true'i2 
■^k ) 
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Figure 8: Kullback-Leibler divergence between the true product of WN densities and the proposed 
approximations. 


between estimates Xk and true state variables We simulated the system for /smax = 100 time 

steps and compared the angular RMSE of all estimators. The results from 100 Monte Carlo runs are 
depicted in Fig. [T^ 

In the scenarios with additive noise, it can be seen that the proposed hlter performs very well 
regardless of the amount of noise. Only the particle filter with 100 particles is able to produce similar 
results. However, it should be noted that the proposed hlter uses just hve samples. The particle 
hlter with 10 particles performs a lot worse and fails completely for small noise as a result of particle 
degeneration issues. Both variants of the UKF perform worse than the proposed hlter. Particularly the 
UKF with two-dimensional state does not work very well, which can be explained by the inaccuracies 
in the conversion of a one-dimensional into a two-dimensional noise noise. 

When non-additive noise is considered, the proposed hlter even signihcantly outperforms the 
particle hlter with 100 particles. As a result of the low number of particles and the associated issues 
regarding particle degeneration, the particle hlter with 10 particles has the worst performance. 

8. Conclusion 

In this paper, we presented a framework for recursive hltering on the circle. The proposed hltering 
algorithms can deal with arbitrary nonlinear system and measurement functions. Furthermore, they 
can be used in conjunction with different circular probability distributions. We have shown that the 
prediction step can be performed based on circular moments only, without ever assuming a particular 
distribution. 

For the purpose of evaluation, we have considered several aspects of the proposed methods. First 
of all, the accuracy of deterministic approximations was evaluated by considering the error when using 
them to propagate a continous distribution through a nonlinear function. We have found that the 
proposed deterministic approximation with five samples yields good results for most practical scenarios. 


22 



















ai = 0.1 


CN 

o 


lO 

d 


S 6 


CT 2 


CO 0.5 


0.04 


S 0.03 


0.02 


a- 0.01 



- VM 



/ ' ' 

1 ' 

/ ' 

/ ' 

V i 


2 4 6 


y 

.'l“ 


l' * 

t ' 


2.5 r 

0) p 
o ^ 
c 
CO 

w 1.5 

■O 

■D 

0 1 
TO 

^ 0.5 
0^ 


cji = 0.4 




0-1 = 1.0 




1 




0.06 


■o 0.04 


• 0.02 


0 0.08 

o 

c 

I 0.06 

'-O 

■p 0.04 

CO 

I 0.02 




M »l 
I U , 

M. 



Figure 9: distance between the true product of WN densities and the proposed approximations. 


Second, we evaluated the novel moment-based WN multiplication method and show that it is superior 
to the previously published method based on fitting a VM distribution. Finally, we evaluated the 
proposed filtering algorithms in several scenarios and compared it to standard approaches. These 
simulations show the advantages of using a circular filtering scheme compared to traditional methods 
intended for the linear case. 

Future work may include extensions of the proposed methods to other manifolds such as the torus 
or the hypersphere. Additionally, consideration of multimodal circular distributions may be of interest, 
for example by means of WN or VM mixtures. 
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